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ABSTRACT 

We have developed an analytical model to describe the evolution of anisotropic galactic out- 
flows. With it, we investigate the impact of varying opening angle on galaxy formation and 
the evolution of the intergalactic medium. We have implemented this model in a Monte Carlo 
algorithm to simulate galaxy formation and outflows in a cosmological context. Using this algo- 
rithm, we have simulated the evolution of a comoving volume of size {12h~^Mpc)^ in the ACDM 
universe. Starting from a Gaussian density field at redshift z = 24, we follow the formation of 
^ 20, 000 galaxies and simulate the galactic outflows produced by these galaxies. When these 
outflows collide with density peaks, ram pressure stripping of the gas inside the peaks may result. 
This occurs in around half the cases and prevents the formation of galaxies. Anisotropic outflows 
follow the path of least resistance and thus travel preferentially into low-density regions, away 
from cosmological structures (filaments and pancakes) where galaxies form. As a result, the num- 
ber of collisions is reduced, leading to the formation of a larger number of galaxies. Anisotropic 
outflows can significantly enrich low-density systems with metals. Conversely, the cross-pollution 
in metals of objects located in a common cosmological structure, like a filament, is significantly 
reduced. Highly anisotropic outflows can travel across cosmological voids and deposit metals in 
other, unrelated cosmological structures. 

Subject headings: cosmology — galaxies: formation — intergalactic medium — methods: ana- 
lytical 



INTRODUCTION 



Galactic outflows play an important role in the evolution of galaxies and the intergalactic medium 
(IGM). Supernova explosions in galaxies create galactic winds, which deposit energy and metal-enriched 
gas into the IGM. These outflows are necessary to explain many observations and to solve many problems 
in galaxy formation, such as the high mass-to-light ratio of dwarf galaxies, the observed metallicity of the 
IGM, the entropy content of the IGM, the abundance of dwarf galaxies in the Local Group, the overcooling 
problem, and the angular momentum problem. 

High-resolution, gasdynamical simulations of explosions in a single object reveal that outflows gen- 
erated by such explosions tend to be highly anisotropic, with the energy and metal-enriched gas being 
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Furthermore, several observations support the existence of anisotropic outflows 
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Indirect support for the existence of a nisotropic outflows comes from the observed enrichment of systems 
around the mean density of the universe (ISchave et al.ll2003t iPieri & Haehnelt"2004') and the enrichment of 
systems far from known galaxies at z ~ 3 (jPieri. Schave. fc Aguirre ,200 6; Songaila 2006). It may be chal- 
lenging to enrich such regions with isotropic outflows even with the inclusion of enrichment from poorly 
understood Population III stars. Early indications are th at Population III stars are unlikely to pollute the 
IGM to a large extent ( Norman. O'Shea. fc PaschosI 2004 ). Anisotropic outflows may also pro vide an expla 



nation for part of the observed scatte r in the metallicity in the IGM, which is still unexplained (jSchave et al 
2003HPieri. Schave. fc Aguirrelbooeh . 



Several simulations of galactic outflows in cosmological contexts have been performed. Typically, these 
simulations use cubic comoving volumes of size ~ (lOMpc)^, containing thousands of galaxies. The methods 
used can be divided into two groups : analytical methods and numerical methods. The analytical meth- 
ods describe the expansion of outflows in simulations that are either Gaussian random realizations of the 
density power spectrum or N-body simulations combin ed with prescriptions for galaxy formation. Out- 



flows are represented using an analytical solution (e.g. IScannapieco fc Broadhurstll200ll hereafter SBOl 



Theuns. Mo. fc Schavel2001 : Bertone. Stoehr. fc White 2005 ) that currently assume isotropy. 



The numerical simulations use hydrodynamical algorithms such as SPH, an d outflows are generated in a 
varie t y of ways: by imparting SPH particles with a large velocity component (e.g. IScannapieco. Thacker. fc Davis 



2001 



Springel fc Hernquist 



SPH particles (e.g. iTheuns et al 



2003 



. .Qppenheimer fc Pavel 200a), depositing additional thermal energy into 
2002b[ ). or taking outp ut from completed S PH simulations and calculating, 



a posteriori^ the propagation of outflows into the IGM (jAguirre et al.ll200ll ). Unlike the analytical methods 
cited above, which assume isotropic outflows, all these numerical approaches have the potential to generate 
anisotropic outflows. However, we believe that there are some limitations to these numerical approaches, 
which motivates us to introduce an analytical model for anisotropic outflows. 

With any SPH simulation, we need to be concerned with the li mited resolution of the algorithm. Consider 
for example the simulations of IScannapieco. Thacker. fc DavisI (|2001l ). These authors identify galaxies of 
radius rjv, and rearrange the SPH particles located between rjv and r^ = 2rM into two uniform, concentric 
spherical shells located at radii 0.9ro and ro, which are then given an outward velocity. The outflow is 
therefore initially isotropic but becomes anisotropic as it propagates into a non-uniform external medium. 
We can see two potential problems with this approach. Firstly, the structure responsible for generating the 
anisotropy might be the galaxy or its environment, in which case the rearranging of particles into concentric 
shells would erase that structure entirely. Secondly, 52 particles per shell (the number they used) provides 
a good covering of the solid angles initially, but as the outflow expands and the particles in the shells move 
apart they become more like individual pressure points pushing on the external medium, and this could lead 
to an artificial mixing of the outflow and the external medium by Raleigh- Taylor instability. 



The approach of lAguirre et al.l (|200ll ) is radically different. It consists of identifying galaxies in an output 
from an SPH simulation and calculating the propagation of outflows from these galaxies in Na different 
directions. Since the resistance encountered by the outflows will be direction-dependent, outflows will start 
isotropic but then become anisotropic as the distance travelled by outflows will vary with direction. There 
are two limitations to this approach. Firstly, since it uses the output of an SPH simulation and introduces the 
outflows a posteriori^ the feedback effect of these outflows cannot be simulated (i.e. outflows do not influence 
the formation of other galaxies). Secondly, in this approach gas elements move radially and encountering 
high-pressure gas will dissipate their energy and rapidly slow down. In the real universe, that gas element 
will likely acquire a tangential velocity component that will redirect it towards regions where the resistance 
of the external medium is weaker. 
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To overcome these various limitations and study anisotropic outflows in a cosmological context, including 
the effect of feedback, we have designed an analytical model for galactic outflows. This can then be combined 
with either an analytical method or a semi-analytic method for the description of galaxy formation. In this 
paper, we use an analytical Monte Carlo method. Calculations performed with an N-body semi-analytic 
approach will be presented in a forthcoming paper (jMartel. Grenon. &: Pierill2007f ). 



This paper is set out as follows. In §2, we describe our analytical model for anisotropic outflows. In §3, 
we describe our Monte Carlo method for cosmological simulations. Results are presented in §4, and their 
implications are discussed in §5. Conclusions are presented in §6. 



2. A MODEL FOR ANISOTROPIC OUTFLOWS 

Consider a density field with a local density maximum at some position, P. A halo may collapse at that 
point and may go on to form a galaxy, which will then produce an outflow. Our goal is to design an analytical 
model for the geometry and evolution of this outflow that takes into account the physical properties of the 
galaxy, the density distribution of the matter surrounding the galaxy, and the global properties of the IGM. 

Using the high-resolution simulations and the observations of galactic outflows as a guide, we will 
represent outflows as "bipolar spherical cones." In a spherical coordinate system (r, 9, 0), the outflow occupies 
the volume deflned by r < R{t), 9 < a/2 or 9 > w — a/2, and < < 27r where R{t) is the radius of the 
outflow and a is the opening angle. This is illustrated in Figure [T] 




Fig. 1. — Geometry of the isotropic and anisotropic outflows. The isotropic outflows are spherical; the 
anisotropic ones are bipolar spherical cones. R{t) is the time-dependent radius of the outflow, a is the 
opening angle, and e the direction of the outflow. 

In the limit a = tt, the outflow becomes spherical, and is entirely described by the radius, R{t). If 
a < TV, the outflow is anisotropic, and two additional parameters must be specified: the opening angle a, and 
the direction of the outflow, which is defined by a unit vector e. Hence, our model of anisotropic outflows is 
a 3-parameter model. 
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One possible way to describe our model phenomenologically is to imagine that the outflow starts 
isotropic, with all the matter moving radially, but the parts of the outflow that encounter high-resistance 
from the external medium acquire a transverse velocity and get redirected into the directions where the 
resistance is weaker. In this context, we see that our model and the model of lAguirre et alj (|2001l) represent 
two limiting cases. In their model, the energy of the gas expanding into dense regions is entirely dissipated. 
In our model, that energy is entirely funneled into the less-dense regions. 



2.1. The Opening Angle 



It is n ot clear what value one sho uld be using for the opening angle. The simulations of lMac Low fc Ferrara 



(1999) and Martel &: Shapirol ( 2001a ) show that the angle actually varies as the expansion proceeds. It starts 



at a low value, a ~ 10° — 45° near the site of the explosion. Once the outflow reaches the low-density (around 
or below the mean density of the Universe) regions it "fans out" , and the opening angle increases to values 
a ~ 45° — 100°, or even larger. The radius also varies with opening angle. Hence, our representation of out- 
flows as bipolar spherical cones of fixed opening angle is a convenient simplification. Until we have a better 
understanding of the morphology of the outflows (something that would require more precise simulations) , 
we will treat the opening angle as a free parameter that can take any value from a = 180° (the isotropic 
limit) to a 20° (the jets limit), keeping in mind that the larger values are likely to provide a more accurate 
description of the outflows once they reach the low-density regions. 

It is reasonable to expect that the opening angle will vary for individual outflows on a case-by-case basis 
due to their differing driving pressure. We treat the chosen value as a typical opening angle with an aim to 
reproducing the global properties of anisotropic outflows for a cosmological volume. 



2.2. The Direction of the Outflow 



The outflows are expected to take a path of least resistance out of a galaxy that forms at the location 
of a density peak. In the current literature, two characteristic scales for this path of lea st resistance have 
been considered: the galactic scale, in which a disk may form ( Mac Low fc Ferrara 1999f) . and larger scale 
filamentary or pancake structures (jMartel fc Shapirol [2001alfbl '). In the former case, outflows are directed 
along the rotation axis of the galaxy. In the latter the pressure of the surrounding medium determines 
the direction of the outflow. We have selected the latter as the defining scale, upon which our anisotropic 
outflows are based, for the following reasons: 



• It is unclear if the galaxy has already formed a well-ordered disk by the time much of the initial 
starburst has occurred. Outflows may start while the galaxy is still in the process of assembling. 
Hydrodynamical simulations of galaxy formations suggest t hat starbursts take place during a chaotic 
merger events at high redshift, before a disk has assembled (jBrook et al.ll2005L 120061 ). 



• Despite recent resul ts suggesting the ir rotation axes tend to lie preferentially parallel to the plane of 
pancake structures ( Noh fc Lee 20061 and references therein) the orientation of these disks are largely 
randomized. Averaging over all outflowing galaxies, only a weak directionality is expected before the 
effect of larger scale structures begins to dominate as the outflows continue to expand. 

• The importance of disk-scale effects is dependent on the locations of SNc within the disk. Bipolar 
outflows due to SNe explosions in a disk are a result of the relatively high pressure gradient along the 
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rotational (minor) axis compared to the major axis. It seems reasonable to assume, therefore, that for 
off-center SNe, where the relative difference in pressure gradient is less acute, the outflow will be less 
bipolar. 



Figure [2] shows an intermediate stage of a simulation ( Martel fc Shapiro! 2001a ) of an explosion inside 
a dwarf galaxy that is forming at the intersection of two emerging filaments inside a cosmological pancake . 
This simulation was performed with an Adaptive SPH algorithm (jShapiro et al.l Il996l : lOwen et al.l Il998[ ) 
with 64^ gas and 64'^ dark-matter particles. The outflow is clearly anisotropic, bipolar, and propagates in 
the direction normal to the pancake plane. Interestingly, the central region where the explosion takes place 
has a nearly isotropic density profile. Hence, it is the anisotropy of the outer regions that results in an 
anisotropic outflow. 

In our Monte Carlo simulations, we assume that galaxies form at the location of local density peaks 
in the matter distribution, and we determine the direction of the outflow by finding the direction of least 
resistance in the vicinity of that peak. We discuss in greater detail the concepts of "peaks" and "vicinity" 
in §§3.2 and 3.3 below. 

Consider a local density maximum located at some point, P, inside the computational volume. We 
center a cartesian coordinate system (x, y, z) on that point, and perform a Taylor expansion of the density 
contrast up to second order, 



5{x,y,z) = 5] 



peak 



Ax^ - By^ - Cz^ - 2Dxy - 2Exz - 2Fyz . 



(1) 



This expression contains no linear terms since the density is a local maximum. In practice, we consider the 
density distribution inside a sphere of radius, R* , centered on the point, P, and perform a least-square fit 
of equation ([1]) to determine numerically the values of the 6 coefficients A, B, C, D, E, F. Once these 
coefficients are determined, we rotate to a new coordinate system (x', y', z'), such that the the cross-terms 
vanish. Equation ^ reduces to 

dix', y', z') = ,5peak - A'x'' - B'y'' - C'z'^ . (2) 



It is actually straightforward to perform this change of coordinates. We combine the coefficients of equa- 
tion IH) to form the following matrix, 

A D E' 

M = D B F . (3) 
EEC 



We then diagonalizc this matrix. Since this matrix is real and symmetric, all three eigenvalues are real. 
These eigenvalues are the coefficients A' , B' , C", respectively, and the 3 corresponding eigenvectors give us 
the directions of the three coordinate axes x' , y', and z', respectively. 

The three coefficients A', B', and C" are always positive, otherwise the point, P, would not be a 
maximum, but rather a saddle point or a minimum. Each coefficient is a measure of how fast the density 
drops as we move away from the peak in the direction corresponding to that coefficient. Therefore, the 
largest coefficient corresponds to the direction along which the density drops the fastest. We will interpret 
this as being the direction of least resistance, and assume that the outflow will naturally follow that direction. 
For instance, if B' happens to be the largest of the 3 coefficients, then the outflow will be directed along the 
y'-axis. 



Fig. 2. — Snapshot of a hydrodynamical simulation of an explosion inside a cosmological pancake. Blue: 
density isosurface, showing the pancake. Notice the circular ripple resulting from the explosion. Red: 
temperature isosurface, showing the outflow. 



2.3. The Expansion of the Outflow 



In this section, we present the equations describing the expansion of the outflow. The technique used 
for solving these equations is shown in Appendix A. 



2.3.1. Basic Equations 



Tegmark, Silk, & Evrard (1993, hereafter TSE) presented a fo rmulation of the expansion of isotropic 



outflows that was based on previous work by numerous authors (ICox fc Smith 



1977; 



Weaver et al 



McCrav fc Kafatos 



1977t 
1987t 



McCrav fc Snow 



Ostriker fc McKeel l 



1974; 



McKee fc Ostriker 



19791: iBruhweiler et al.lll980l: iTomisata. Habe. fc Ikeuchilll980l: 



19881 ) and refined it to account for the expansion of the universe. 



In this formulation, the injection of thermal energy produces an outflow of radius, R, which consists of a 
dense shell of thickness R6 containing a cavity. A fraction, 1 — /„i, of the gas is piled up in the shell, while 
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a fraction, /„, of the gas is distributed inside the cavity. We normally assume (5 ^ 1, ^ 1, that is, most 
of the gas is located inside a thin shell. This is called the thin-shell approximation. 

T he original equations of TSE have b een refined by several authors ( Madau. Ferrara. fc ReesI 2001 



SBOl; Scannapieco. Ferrara. fc Madaull20o5 hereafter SFM) to include various additional physical processes 



Hence, there is not a "unique" form of the equations. In this paper, the evolution of the shell radius, R, 
expanding out of a halo of mass, M, is described by the following system of equations, 

L 5Rp 

(5) 



27ri?3 R ' 

where a dot represents a time derivative, fi, fib, and H are the total density parameter, baryon density 
parameter, and Hubble parameter at time, i, respectively, L is the luminosity (discussed in §2.3.3), p is the 
bubble pressure resulting from this luminosity, and Pext is the external pressure of the IGM. The four terms 
in equation (|4|) represent, from left to right, the driving pressure of the outflow, drag due to sweeping up the 
IGM and accelerating it from velocity HR to velocity R, and the gravitational deceleration caused by both 
the expanding shell and the halo itself. The two terms in equation ([5]) represent the increase in pressure 
caused by injection of thermal energy, and th e drop in pressure caused by the expansion of the outflow, 
respectively. For details, we refer the reader to Ostriker fc McKee ( IQSSh and TSE. 



These expressions most closely resemble that of SFM, the one difference coming fr om their use of a radi 



ally d ependent mass in the last term of equation (j4|) , derived by assuming a NEW profile (jNavarro. Frenk. fc White 



19971 ). This is motivated by the gravitational drag on a spherical outflow expanding out of a spherical col- 
lapsed halo, and is no longer an adequate description in the case being examined since an anisotropic outflow 
will more rapidly escape an elliptical halo. The influence of the gravitational potential well will be an inter- 
mediate case between this and a point mass. In practice we take the simpler form by assuming a point mass, 
as in SBOl. We justify this approximation as follows: early in the outflow lifetime the effect of gravity on the 
outflow is negligible compared to other effects (see Appendix A). As the outflow expands, the gravity terms 
eventually become important, by which time most outflows have left the halo from which they originate. 
Notice that the original formulation of TSE did not include the term —GM/R'^, and also ignored the external 
pressure Poxt- 

These equations provide a full description of isotropic outflows. We now modify these equations to 
describe anisotropic outflows with any value of the opening angle, a. We impose the condition that the 
modified equations reduce to equations ([4]) and ([5]) in the isotropic limit a — 180°. 

First consider equation Q . The two terms in the right hand side correspond to the increase in pressure 
caused by thermal energy injection, and the decrease in pressure caused by the expansion of the outflow, 
respectively. At this point, we need to backtrack. Equation ([5]) was derived from the following system of 
equations (TSE): 

Et = L-pV, (7) 



where Et is the thermal energy inside the outflow, and V is the volume of the outflow. We take the time 
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derivative of equation ([6]), and eliminate Et using equation We get 

. _ 2L 5pV 



(8) 



If we now substitute V — AttR^/S, we recover equation ([5]). However, a bipolar anisotropic outflow with 
opening angle, a, has a volume given by 

F=^(l-cos-). (9) 

We substitute equation ([9| in equation dH), and get 

L 5Rp 



27ri?3[l-cos(a/2)] R 



(10) 



Physically, as a decreases, the energy is injected in a smaller volume, resulting in a larger pressure, p. 
This equation replaces equation ([5|). In our model, this is the only modification we need to make to the 
equations describing the expansion of the outflow (although a modification to the luminosity will be required). 
Equation ^ remains unchanged. 



2.3.2. The External Pressure 



The external pressure, Pcxt, depends upon the density and temperature of the IGM, which can be quite 
complex. But in the context of the a nalytic approximation we are us ing, we are justified in making some 
simplifying approximations. Following iMadau. Ferrara. fc ReesI (j200ll) . we will assume a photoheated IGM 
with a fixed tempe rature, Tigm = lO'^K . We will also assume an IGM density, pigm, equal to the mean 



baryon density, pb. iHui fc GnedinI (119971 ) performed detailed hydrodynamical calculations of the equation of 
state for a photoheated IGM, and showed that the temperature is a function of the density. But according 
to their Figure 1, a temperature of lO^K is appropriate for de nsities, p ~ p^,. It should also be noted that 



the temperature of the IGM is redshift-dependent as shown by ISchave et al.l (j200Cll ) and we will neglect this 
here. Our choice of fixed external temperature and density will be discussed further in §5. The external 
pressure is then given by 

PbkTiGM 30h^o^fofc7iGM(l + z)^ 



PeKtiz) 



SnGp 



(11) 



where z is the redshift, p is the mean molecular mass, and subscripts indicate prese nt values. The val ue of 
p depends upon the ionization state of the gas. We assu me that hydrogen is i onized (Page et al."2006') and 
helium is singly-ionized in the redshift range 3 < z < 7 (jTheuns et al.ll2002at IChoudhury fc F errara 2005), 
which is most of the redshift range in which our outflows are actively expanding. Also, even when much of 
the IGM is still neutral, it is reasonable to expect the source galaxy to have photoionized the region ahead 
of the outflow. We will assume that this is the case throughout the simulations. The mean molecular weight 
is the n p — 2/{4 — 3Y) a.m.u. (atomic mass units). For a helium abundance, Y = 0.242 ( Izotov fc Thuan 
20041 ). we get p = 0.611 a.m.u.. 
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2.3.3. The Luminosity 

The luminosity, L, is the rate of energy deposition or dissipation within the outflow and is a sum of five 
terms, 

L{t) — LsN — icomp — ino2 " -^ion + ^diss , (12) 

where Lsn is the total luminosity of the supernovae responsible for generating the outflow, Lcomp represents 
the cooling due to Compton drag against CMB photons, Lno^ represents the cooling due to 2-body interac- 
tions. Lion represents the cooling due to ionization, and Ldiss represents heating from collisions between the 
shell and the IGM. Following SFM, we will assume that the first two terms, supernovae and Compton drag, 
dominate, and will neglect the remaining terms. The supernovae luminosity, for a galaxy forming in a halo 
of mass, M, is given by 

tburst Afrcq ^^^* ^ JIq ^ ^ IM© ^ ® ' ^ ^ 

where /w is the fraction of the total energy released that goes into the outflow, /, is the star formation 
efficiency, and the total mass in stars formed during the starburst is M* = /,Af ilb^o/^^o- Afreq is the mass of 
stars required to fo rm one SN and we take a value of Afj-cq — 89.7M0 (derived using a broken power-law IMF 
from Kroupall2001[ ) and the energy released by each of these SNe is — 10^^ ergs. iLeitherer fc Heckman 
(jl995r ) show that for instantaneous star formation, as considered here, the SN outflow phase associated with 
this burst last for a duration of tburst = 5 x 10'' yr and the delay before the first SNe is negligible. We use 
the galaxy mass dependent determination of /w given by SFM, /w(Af) — Q ."ib b{M) j E b{M = 2 x IO^Mq), 
where 

r 1.0, Nt < 1; 

Sb{M) = I 1.0 + 0.165\n{Nt-^), 1 < < 100; (14) 

[ [1.0 - 0.165 In(lOO)] 100iVt"\ 100 < Nt] 

and Nt = 1-7 x 10~'^{nbfi/flo)M/{lMQ). 

The ionization state inside the outflow is even more uncertain and may be a function of distance from 
the source as the gas may recombine as it expands from its source galaxy. Until we have a more precise model 
for the internal structure of the outflow, we will simply make the same assumption as for the surrounding 
IGM: ionized hydrogen and singly-ionized helium. However, it is worth noting that the extreme case, where 
outflows contain fully ionized gas, involves only a minor modiflcation of the Compton drag. Following TSE, 
the Compton luminosity is given by 

15 ymeC^ J [he) 



Lcomp = — {(Ttcne) ( _ ) {kT^)'^j^-^ , (15) 



where at is the Thompson cross section, is the CMB temperature, and Ue and Te are the electron number 
density and temperature inside the bubble, respectively. The pressure, p, inside the bubble is given by 

p = nkTe = {ne + nion)kTeV . (16) 

We assume that hydrogen is ionized and helium is singly- ionized and hence riion = ne, which is also the 
approximation made in TSE. Equation (|16p becomes UekTe = p/2. We now eliminate n^kTe in equation (|15p . 
We also eliminate V using equation ©, and get 

^---!^S(1f)'^"^f)<'-''''«'' '''' 



- 10 - 



where T^o is the present CMB temperature. We used — T^q{1 + z), which is vahd over the range of 
redshifts we consider. 



2.3.4- The Final Stage of the Outflow 

The expansion of the outflow is initiaUy driven by the supernovae luminosity. After the supernovae 
turn off, the outflow enters the "post-supernova phase." The pressure inside the outflow keeps driving the 
expansion, but this pressure drops since there is no energy input from supernovae. EventuaUy, the pressure 
wiU drop down to the level of the external IGM pressure. At that point, the expansion of the outflow will 
simply follow Hubble expansion. Since we assume identical mean molecular weights inside and outside the 
outflow, the condition p = pcxt implies T = Tjgm = lO^K. 

Figure [3] shows the isotropic outflow solution for a halo of mass, M = 2.8 x lO^M© that collapses at 
redshift, z = 8.10, in a ACDM universe. The cooling time is short, 3Myr, so the galaxy forms soon after, at 
z = 8.07. The outflow turns on immediately, as discussed in the previous section, and goes on to form the 
largest outflow in our simulations. 

The top panel shows the total luminosity, and the middle panel shows the comoving radius. The bottom 
panel shows the pressure, p, of the outflow and the external pressure, Poxt, of the IGM. Initially, during the 
supernova phase, the outflow is driven by the energy input from supernovae and the radius steadily increases. 
During this phase, the supernova luminosity, Lsn, exceeds the Compton luminosity by factors of several and 
for this particular outflow dominates it entirely. The pressure diverges at the redshift of supernova turn-on 
(because we neglect the flnite volume of the region containing the supernovae), but then drops as t~^^^ 
(see Appendix A). At z = 7.63, the supernovae burn out, and the only contribution to the total luminosity 
is the energy loss by Compton drag. The outflow enters the post-supernova phase. The turn-off of the 
supernovae produces a change of slope in the pressure, which in turn produces a change of slope in the 
Compton luminosity. The post-supernova phase lasts until redshift z — 4.95. At that point, the pressure of 
the outflow becomes equal to the external pressure, and the outflow switches to pure Hubble expansion. 



3. THE MONTE CARLO METHOD 

We use essentially the same Monte Carlo method as in SBOl, generalized to include anisotropic outflows. 
In our dynamical equation for the evolution of the outflow [eq. (|4])], we include the external pressure, Pext, 
an effect that was neglected in SBOl, but included in SFM. 



3.1. Cosmological Model and Initial Conditions 

We consider a ACDM model with present density parameter, Hq = 0.27, baryon density parameter, 
ilbfl — 0.04444, cosmological constant, Aq — 0.73, Hubble constant, Hq = Tlkms^^Mpc^^ {h = 0.71), 
primordial tilt, n.s = 0.93, and CMB temperature, Tcmb — 2.725, consistent with the results of WMAP 
(jBennett et al.l 120031 ). We simulate structure formation inside a comoving cubic volume of size, Lbox = 
12/i~^Mpc, with periodic boundary conditions. To generate initial conditions, we lay down a cubic grid of 
size 512 X 512 x 512 in the computational volume, and compute the initial density contrast, Si, at initial 
redshift, Zi ~ 24. We also compute, on 10 similar grids, the same density contrast, but filtered at the 10 



- 11 - 



7.2 



7.15 



^ 7.1 
o 



7.05 
^ 0.3 



^ 0.2 
c 

o 0.1 

o 
o 

K 



O 



o 



12 



14 



-16 




H h 



post-SN 



H ^ h 



H ^ h 



Hubble 



H ^ h 




Fig. 3. — Evolution of the largest outflow. The middle panel shows the comoving radius of the outflow versus 
redshift, z. The thick lines separate the various phases of the expansion (supernova-driven, post-supernova, 
Hubble expansion). The top panel shows the total luminosity, L = Lg^ ~ Lcomp, versus redshift. The bottom 
panel shows the external IGM pressure (solid line) and the pressure of the outflow (dotted line) . 



mass scales: Mi, M2, . . . , Miq. As we use a Gaussian filter, the mass, M, and comoving radius, Rf, of the 
filter are related hy M — {27r)^/^R^po, where po — S^oHq/SttG is the present mean density of the universe. 
The values of the mass and length filtering scales are listed in Table [TJ The last two columns indicate the 
filtering radius in units of the grid spacing. A, and in units of the box size, Lbox- 



The method for generating these 11 grids is described in great detail bv iMartell (j2005l ). Essentially, we 
work in /c-space, by generating, on a 512'^ grid, the density harmonics, S{k), corresponding to a ACDM power 
spectrum at z = Zi — 24. Once the density harmonics are generated we take the inverse Fourier transform 
to obtain the initial density contrast. Si. This gives us our first grid in real space, with the unfiltered density 
contrast. To get the 10 filtered grids, we first multiply the density harmonics by the Fourier transform of 
the filter, and then take the inverse Fourier transform to get the filtered density contrast. Notice that in 
our simulations the value of the initial redshift, z^, is actually irrelevant, as long as it is larger than the 
collapse redshifts of all the density peaks that are resolved in our simulations. It turns out that the first 
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peak collapses at redshift, z = 18.55. 



3.2. The Local Density Peaks 

For each of the 10 filtered grids, we identify the location of the local density peaks. These peaks are 
defined as grid points where the density contrast is positive, and its value exceeds the values at the 26 
neighboring grid points (taking periodic boundary conditions into account). For each density peak on each 
filtered grid, we compute a collapse redshift, Zcoii, and a direction of least resistance, e. The collapse redshift 
is obtained by solving numerically the following equation, 

<5c(zcoii) = ^.%%2^, (18) 



where 5i is the initial density contrast at the peak, and ^+ is the linear growing mode (for flat, A 7^ 
models, see, e.g.. lMartel 199L equation [18]). The critical value, (5c, is equal to 1.69 for an Einstein-de Sitter 
model, and slowly varies with redshift for other models. In our simulations, we simply assume be — 1.69. If 
equation (|18p has no solution, we are simply dealing with a density peak that would collapse in the future, 
and we just ignore it. 

To check for consistency, we performed a N-body simulation of structure formation in a box of the same 
size, for the same cosrnologi cal model and density fluctuation power spectrum. We used a P'^M algorithm 
([Hocknev &: Eastwoodlll981r ) with 256'^ equal-mass particles. The mass per particles was 1.087 x 10 ''M0, 



hence the lowest mass scale M\ corresponds to 7 particles which, though insufhcient for determining the halo 
density pr ofile, is fine for th e simple locating of halos that follows. We used a standard friends-of-friends 



algorithm (jPavis et al.lll985l ) to identify halos at various redshifts between z = 24 and z = 2. For this, 
we used 2 different linking lengths. We first use the "standard" value I — 0.25Aa::, where Ax is the mean 
particle spacing. We also used a value I — (187r2)~i/'^Ax — 0.1779Aa;, corresponding to a density increase 
by a factor of Vi'n^ . This is consistent with the assumption that collapsed, virialized halos have a density 
equal to 187r2 times the mean background density (see eq. [19] below). Only halos containing 6 particles or 
more are included. This corresponds to a minimmn mass of 6.522 x IO^'Mq. 

In Figure 21 we plot the number of collapsed halos present in the N-body simulation versus redshift 
(dashed lines). We find more halos with / — 0.25 than with / = 0.1779. This might seem surprizing since 



Table 1. Mass and Length Filtering Scales 



Filter name M [Mq] i?/ [kpc] i?//A i?//Lb, 



MOl 


7.61 


X 


10^ 


50.4 


1.53 


0.00299 


M02 


2.53 


X 


108 


75.2 


2.28 


0.00445 


M03 


8.42 


X 


108 


112.3 


3.40 


0.00664 


M04 


2.80 


X 


109 


167.6 


5.08 


0.00992 


M05 


9.32 


X 


109 


250.2 


7.58 


0.0148 


M06 


3.10 


X 


lOlo 


373.6 


11.31 


0.0221 


M07 


1.03 


X 


10" 


557.8 


16.90 


0.0330 


M08 


3.43 


X 


10" 


832.7 


25.23 


0.0493 


M09 


1.14 


X 


1012 


1243.1 


37.66 


0.0736 


MIO 


3.80 


X 


10^2 


1855.9 


56.22 


0.110 



- 13 - 




Fig. 4. — Number of collapsed halos present in the computational volume, versus redshift. Long and short 
dashes: N-body simulation with friends-of-friends algorithm, with linking lengths, I — 0.25Aa; and / = 
0.1779Ax, respectively; dotted line: Press-Schechter approximation; so^id line: Monte Carlo simulations. 



halos "break up" as I is reduced. However, many halos fall below the threshold of 6 particles when I is 
reduced. The two methods differ by about 20% at low redshift, but the difference increas es at high redshift. 
The dotted line shows a calculation based on the Press-Schechter (PS) approximation (jPress fc Schechter 



1974I ). The agreement with the N-body simulations is excellent at low redshift, while at very high redshifts 



{z > 15) the PS approximation overestimates the number of halos. The solid line shows the number of 
collapsed halos in our Monte Carlo simulations. This was obtained by counting, at given redshifts z, the 
number of halos with Zcoii > and removing the halos that have been destroyed by mergers (see §3.4 
below). The number of halos is lower for the Monte Carlo simulations than for the N-body simulation 
or the PS approximation, at almost all redshifts. We interpret this result as follows: in the Monte Carlo 
simulations, only the matter located in overdense regions can eventually end up inside collapsed halos, while 
in the N-body simulations, all the matter in the system can eventually end up in halos. Notice also that 
the comparison between the PS and N-body results, and the Monte Carlo results is complicated by the fact 
that the mass spectrum of halos is discrete for the Monte Carlo simulations. For the N-body simulation and 
the PS approximation, the choice of an appropriate minimum mass is not obvious, and the results are quite 
sensitive to that choice, especially at high redshift when most halos have a small mass. 
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These results also vindicate our choice of an initial redshift of, Zi — 24. Clearly the system is still in 
the linear regime at that redshift. The very first halo collapses at redshift z = 18.55 in the Monte Carlo 
simulations, and at rcdshifts, z ~ 16 — 17, in the N-body simulation, while the PS approximation predicts 
that a computational volume of that size will typically contain one halo at z — 19. Indeed, the probability 
of finding one halo of mass, M = Mi = 7.61 x lO^M©, at z = 24, according to the PS approximation, is 
about 1/1000. 

3.3. The Direction of Least Resistance 

We compute, for each peak, the direction of least resistance, using the method described in ^2.2^ The 
radius, i?*, defines the region around the peak over which the least-square fit to equation ([T]) is performed. 
Its value will affect the determination of the coefficients. A, B, C, D, E and F, and ultimately the direction, 
e, of the outflow. The largest possible value of R* should be the filtering scale, since this scale corresponds 
to the physical extent of the peak (using a larger value would amount to including matter that belongs to 
different peaks). As for the minimum value, R* must clearly be at least equal to the grid spacing, A. It 
turns out that this condition is insufficient. If R* < 2^/'^A, only the 6 nearest grid points to the peak are 
included in the fit. In that case, the matrix, M, is already diagonal, and the only directions allowed for the 
outflow are along the unrotated coordinate axes x, y, z. We must at least include the 12 next nearest grid 
points, located at a distance equal to from the peak. 

To test how the direction of the outflow depends on the particular choice of R*, we considered the filters, 
M03, M04, . . ., MIO, and computed, for all peaks, two values of e: one using R* = 2A, and one using R* 
equal to the filtering scale, Rf. Figure [5] shows histograms of the angle in degrees between these two values 
of e. We excluded from this test the mass filters MOl and M02, because the range of allowed values for R* 
becomes very narrow, and there is essentially no uncertainty in these cases. 

For the vast majority of peaks, the angle is below 4 degrees and it exceeds 10 degrees only for very few 
peaks. The interpretation of these results is that our method for determining the direction of the outflow 
has an uncertainty of a few degrees, resulting from the ambiguity in the choice of R* . In the limit when the 
opening angle of the outflow is much larger than that uncertainty, the region of space containing the outflow 
is essentially unchanged. We will consider opening angles of several tens of degrees, while the uncertainty in 
the direction of least resistance is of order a few degrees. The final results might vary by a few percent as a 
result, but this is not sufficiently large to be the dominant source of error. 

From this test, we conclude that our technique for determining the direction of least resistance is robust, 
and consequently the particular value chosen for R* is not important. Hence, we decided to set R* equal to 
the filtering scale, Rf. 

Once this is all done, we have, for each of the 10 filtering scales, a list of local density peaks, and for each 
of these peaks, a position, r, a mass, M (the mass of the corresponding filtering scale), a collapse redshift, 
Zcoiij and a direction of least resistance, e. In the model, each peak corresponds to a density fluctuation of 
mass, M (the flltering mass), and radius, Rf (the flltering radius). If left alone, this peak will collapse at 
redshift, Zcoii, to form a halo of mass, M, located at position, r. This protogalaxy will eventually turn into 
a galaxy (by forming stars), and generate an outflow propagating along the direction of least resistance, e. 
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Fig. 5. — Histograms of the fluctuations in the direction of the outflows resulting from variations in R*. 
Notice that the last two bins correspond to the ranges [10,15] degrees and [15,90] degrees, respectively. The 
labels in the panels indicate the mass filtering scale and the number of peaks. 



3.4. Formation of First-Generation Galaxies 



The 10 filtered grids represent the same region of the universe. Hence, we cannot simply assume that 
every peak collapses to form a protogalaxy, since this would violate conservation of mass. We must deal 
with the issue of halos inside halos|^ Consider two density peaks a and b of masses, Ma < Mb (these 
peaks are therefore on different grids). If the separation, |ra — r;,|, between the peaks is smaller than the 
filtering scale, of the larger mass, the two halos that these peaks will form cannot coexist, since the 

smaller halo is inside the larger one. There are then two possibilities, depending on the collapse redshifts. If 
{zcoiijb > (-Zcoii)oj then halo a will never form, since all the matter it contains will be incorporated into halo 
b, that forms earlier. While this does occur it is not common, since in a CDM universe smaller peaks tend 
to collapse earlier than bigger ones. The second possibility is (zcoiOa > {zcoi\)b- In this case, halo a will form 



Throughout this paper, we shall use the term "halo" to designate a collapsed object, whether it is a galaxy or a protogalaxy. 
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first, at redshift (zcoiOa, but will be destroyed later by a merger event at redshift, (zcoiOb, that results in the 
formation of halo b. So, even though halos remain at fixed locations in this simplified Monte Carlo model, 
non-dynamic merger events are included. 

The collapse of a peak leads to the formation of a halo, which we identify as a protogalaxy. The gas 
inside that halo is at the virial temperature, given by 



0.009 K / 6.8 



13 \hX + i 



Mh 



2/3 



(1 + Zcou) 



17o 


1/3 


Ac(Zcou) 


1/3 


_VL(Zco\l)_ 




187r2 





(19) 



(|Eke. Cole, fc Frenklll996l ). where X = 0.76 is the hydrogen mass fraction, and is the ratio of the halo's 
mean density to the critical density. Since this temperature is usually much too high to allow the formation 
of stars, the gas must cool down to a temperature, T <C Tyjr, before s t ars fo rm and the protogalaxy becomes 
a galaxy. As in SBOl, we use the cooling model of I White fc FrenkI (|l991l ). This model assumes that the 
cooling proceeds inside-out. The gas inside a "cooling radius," TcooI, has cooled, while the gas outside rcooi 
remains hot. As TcooI increases with time, the mass, Mcooh of cool gas increases, until all the gas has cooled. 
The mass, Mcooi, evolves according to 



dM, 



cool 



dt 
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3/2 
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(20) 



where A23 is the cooling rate in units of 10~ ergs s~ c m~'^, which i s a fun ction of temperature and metallicity, 
Z, and /o ~ 0.8 [a correction factor suggested by Somervil"l3 ( 1997 ) to bring this result into line with 
more complex analyses]. We determine the required cooling rate using MAPPINGS 11%, a successor to 
MAPPINGS II outhned in lSutherland fc Dopital | 19931 ). We integrate both sides between t = Q and t = icoob 
the cooling time. For a halo of mass, M, and gas mass, Mgas = ^bfiM /^q, we get 



^cool — 



1 



576/o^ VIM. 



M 



^) [A23(r.ir,Z)]-Vr. 



(21) 



Using the solution of the Friedman equation for our ACDM model, we compute, for each halo, the 
epoch, tco\i, at which collapse occurs from the collapse redshift, Zcou. We then add the cooling time to 
get the epoch at which galaxy formation occurs, t^t = tcoW + icooi, and compute the corresponding galaxy 
formation redshift, Zgf. The galaxy formation epoch may be later than the current epoch, in which case 
stars will not form and the halo will remain a protogalaxy, unless metal deposition from outflows modifies 
the cooling rate (see below). 

We refer to these galaxies as "first-generation" galaxies since at the time of formation they are untouched 
by impinging outflows from other galaxies. For these galaxies we assume a cooling rate based on zero 
metallicity. In the following section we will discuss galaxies for which this is not the case. 

We are now ready to proceed with the simulation. Starting at initial redshift Zi — 20 (since there are 
no collapsed halos at higher redshift), we evolve the system forward by redshift steps of Az = —0.005. As 
a result we take 3600 steps to reach the end point for our simulations, which is z = 2. This end point is 
chosen for comparison with Quasar absorption spectra data from z = 2 — 6. Also, our box would no longer 
represent a cosmological volume at lower redshifts. 



^http: / /www. mso.anu.cdu.au/~ralph/map. html 
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At the location of peak, a, a halo forms at redshift, (zcoiOa; unless that peak was located inside a larger 
peak that collapsed first, as described above. During the redshift interval, [(zcoiOa, (^gf)a], while the gas is 
cooling, the halo might be destroyed by a merger. This will happen if there is a peak b for which 

Mfc > , (22) 

(^colOa > (Zcoll)fc > {Zgf)a , (23) 

\ra-rb\<RfM- (24) 

In this case, the halo was able to form, but did not manage to cool before being destroyed, hence the actual 
galaxy never formed. 

If the galaxy does manage to form, we neglect the formation time and life span of massive stars and so 
the newly formed galaxy immediately produces an outflow centered at the halo position, , and propagating 
along the direction, Gq. We evolve this outflow with time, using the solutions given in Appendix A (with 
timesteps corresponding to our redshift steps, Az = —0.005). 

Left alone, the outflow will go through successive phases of driven expansion by supernovae and near- 
adiabatic expansion driven by pressure alone before joining the Hubble flow. However, during either of these 
stages, the galaxy producing the outflow might be destroyed by a merger. In all such cases the outflow 
reverts to a Hubble expansion phase. 



3.5. Subsequent Galaxy Formation 

After the first outfiows are formed, some will begin to strike other density peak^ and modify the way 
these systems evolve. They may enrich these systems with new metals, possibly modifying their cooling time. 
They may also expel the gas content of the object by ram pressure s tripping or shock-heating of the gas. 
Of th ese two mechanisms, the former (stripping) is the dominant one (jScannapieco. Ferrara. fc Broadhurst 



20001 ) and so we will neglect the latter. 



3.5.1. Stripping 

Density peaks may be stripped of their baryons when the swept-up shell of intergalactic and interstellar 
gas incident upon them imparts sufficient momentum on the baryons such that they escape the potential 
well of the peak, i.e. 

^ MoVo > MfcWosc, (25) 

where the mass of the swept up shell is Mq at radius, R, (as in TSE), Vo is the outfiow velocity. Mi, is the 
baryonic mass of the density peak struck, and Vcsc is the escape velocity for successfully stripped gas. I 
is the comoving radius of the collapsing density peak. Using the spherical approximation for gravitational 
collapse, this scale is determined from the non-linear density contrast, 5nl using I = i?/(l + 5NL)"^/^ where 



^ In the following text we will use the term "peaks" to refer to density peaks that will go on to form halos by 2 = 0. 
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1 + (^NL = 9(0 — sin6')^/[2(l — cos^)^]. The parameter, 9, is in turn given by 

0+(Zco\l) 



-sin 61 5+(z) 
"2^ 



where S+{zcoii) and (^+(2;) are the hnear growing modes at collapse and when the outflow strikes the peak 
respectively. 

We assume that stripping may only occur for systems that have not collapsed since the cross section 
to impact is small and would result i n negligible momentum transfer. This approach is supported by 
Sigward. Ferrara. fc Scannapiecol ( 20051 ) who investigate this issue with a numerical analysis of an individual 



object struck by a shock. It is also worth noting that we deal with baryonic stripping in the same way 
whether or not the source galaxy is within the filtering radius of the density peak struck (although it is less 
likely to occur due to the larger escape velocity and mass of the halo struck). 

When the criterion for stripping is met, the density peak is rendered free of baryons and, while it may 
collapse, it will not form a galaxy. If the shell of the outflow does not strip the peak then the peak is enriched 
by the metal-rich gas, which fills the outfiow. 



3.5.2. Metal Enrichment 



Metals are propagated throughout our simulations within the hot bubbles of outflows. Peaks that have 
not been struck by outflows are assumed to have metallicity of [Fe/H] = —3, which is negligible for the 
purpose s of calculating the halo cooling time. Once galaxies are formed metals are produced at rate of 2Mo 
per SN iNagataki fc Satol l Il998r ). Hence the mass of metals in the outflow is 



Mz = U. 



(27) 



where /esc is the fraction of I SM gas blown out with the out flow. We use the value /esc — 0.5 taken from 



Mori. Ferrara. fc Madaul \200i ). This mass of metals is distributed evenly 



the numerical simulations of 
throughout the volume of the outflow. 

It is notable that, in this model, the ratio of the mass lost rate to star formation rate is approximately 
c//* = 5. This is consistent with observations that indicate this ratio is of order 1 for local starburst 



galaxies (e.g. Martin 1999 and Heckman et al. 2000[) . infrared-luminous galaxies (e.g. Rupke et al. II2OO5,) 



and high- 2: galaxies (e.g. 



Pettini et al 



200' 



J" 



When an outflow strikes an uncoUapsed density peak and does not strip it, it modifies its metal content. 
It deposits a fraction of its metals, /depVi^veriap/Kjutflow, where /dcp is a mass deposition efficiency, which we 
set at /dop = 0.9, Vovcriap is the volume of overlap of the uncoUapsed density peak of radius I, and Voutflow is 
the volume of the outflow. This volume of overlap is calculated based on various geometric approximations 
to the overlap of two spherical cones and a sphere, depending on their relative size and orientation. Once 
the halo has collapsed, the addition of metals results in an increase in cooling rate and so a reduction in the 
cooling time, expediting galaxy formation. 

We assume that the SNe in our simulations are only of the Type II variety since these SNe explode 
together shortly after the starburst and lead to a coherent extragalactic o utfiow that will reach large dis- 
tances. As a result we use an alpha-element-rich yield of metals [quoted in Sutherland fc Dopita ( 19931 ) as 
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"primordial" abundance ratios] and distribute our mass of metals accordingly. This provides a value for 
[Fe/H] used to calculate the new cooling rate. 

There is a special case to consider where a source galaxy is embedded in the density peak that its 
outflow strikes. In this case we assume the volume overlap of the outflow is 100% and so it will deposit a 
fraction, /dop, of its metals onto this structure. The main motivating factor behind the inclusion of this mass 
deposition efficiency is to avoid the implausible scenario that an outflow may continue to grow having left 
all its metals behind. The value it should take is unknown; however, the simulations are not sensitive to this 
factor, as we point out below. 

4. RESULTS 

4.1. The Impact of Varying Opening Angle 

We have constructed a model for the evolution of anisotropic outflows and a test bed with which to 
investigate its significance and we are now ready to discuss what this approach tells us. 

Taking an extreme case of anisotropy as an illustration, consider an opening angle of only 40°. At the 
natural end point of our simulations, z = 2, Figure [6] shows a slice of thickness 0.4/i^^Mpc through our 
simulation box. The position and extent of our galactic outflows are indicated along with the position and 
pre-coUapse radius of peaks (i.e. the filtering scale, Rf, of the peak) that will collapse by z = 0. These 
outflows extend a significant distance from their source galaxies and often strike other peaks. Where peaks 
are arranged in a row in the plane of the slice, the outflows appear to favor a direction perpendicular to 
this structure (see zoom in). This is to be expected: the locations of these galaxies trace the underlying 
structure of the dense filament out of which they form by fragmentation, and the outflows follow the path 
of least resistance away from that filament. 

Figure [7| shows the counts of various quantities in units of (/i^^Mpc)^^ by 2; = 2, for varying opening 
angle. The number density of peaks that are hit by outfiows that will go on to form halos by 2; = are 
shown. This decreases as approximately a power law for increasingly anisotropic outflows, which can be 
explained by two factors. Anisotropic outflows take a path of least resistance out of highly overdense regions 
and so will encounter fewer regions that go on to form halos. Also the total volume occupied by outflows is 
not conserved for varying opening angle since anisotropic outflows tend to overlap less (see above) and the 
volume per outflow is not constant. The degeneracy between these two effects is resolved in Figure [5] and 
the accompanying text. The number density of peaks stripped of their baryons is also shown and indicates 
that the fraction of peaks hit that are then stripped is roughly constant (^ 50%). 

The number of galaxies formed by z = 2 rises smoothly for decreasing opening angle in Figure [7] as a 
result of the fall in stripping. We find that metal enrichment of peaks has a negligible impact on the number 
of galaxies; however, those galaxies that form are born slightly earlier. Not all peaks that are hit by outflows 
but not stripped will go on to produce galaxies by z = 2. This is because, despite being enriched by metals 
and collapsing, their baryons do not always cool to form stars by this redshift. As a result, the incidence of 
stripping decreases faster than the number of galaxies increases. This is demonstrated with the sum of these 
two number densities in Figure [3 which decreases with increasingly anisotropic outflows. 
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Fig. 6. — A slice through the simulation box of thickness 0.4ft,~^Mpc at z = 2. Density peaks are shown as 
filled circles of diameter corresponding to their extent prior to collapse. The outflows have an opening angle 
of 40°, and their location and physical coverage are indicated as red wedges. Two zoom-in regions of size 
1.5ft.~"'^Mpc show regions of interest, where galaxies that formed out of a common filament produce outflows 
that are aligned. 
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Fig. 7. — The number density of various objects in units (ft.^^Mpc)"'^ at z — 2, versus opening angle. Solid 
line: number of peaks hit that will go on to collapse before z — 0; dotted line: number of those peaks that are 
then stripped of their baryons; short dashed line: number of galaxies; long dashed line: sum of the number 
of galaxies and stripped peaks. 

In Figure [8] we show various statistics at z = 2 in the simulations: the average distance, i?strip, traveled 
by outflows before stripping occurs, the maximum distance, i?max, traveled by an outflow, the estimated 
volume filling factor of outflows, and the ratio of the number density of hits to the volume filling factor in 
units (h^^Mpc)^^ . The mean distance traveled by outflows before stripping occurs is fiat until the opening 
angle as low as 100° where it begins to increase. As the anisotropic outflows travel farther with decreasing 
opening angle, they fail to hit new peaks as they have already escaped their highly overdense environments. 
It is only when a — 100° that they begin to reach the next dense structure. 

In the top right panel we show the radius of the largest outflow for each opening angle. We also show 
a curve for constant volume using equation ([9]), which gives us i? oc 1/[1 — cos (a/2)]^/'^. This shows that 
while the radii of outflows grow with decreasing opening angle they do not grow fast enough to conserve the 
volume per outflow. 

The estimate of the volume filling factor is calculated based on the sum of the volume of all outflows. 
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Fig. 8. — Top left: mean radius of outflows when they strip peaks of their baryons; top right: the solid line 
shows the radius of the largest outflow and, for comparison, the dotted line shows i? oc — cos a/2] 
which is the relation expected if the volume enriched were conserved for varying opening angle. The estimated 
volume filling factor is shown in the bottom left panel; bottom right: the ratio of the number density of peaks 
hit to the volume flUing factor in units of (h^^Mpc)^^ . 



with a correction factor for the increasing probability of overlap for large values, 



F — 1 — exp 



(12/i-iMpc)3 



(28) 



where Vi is the volume per outflow given by equation © and No is the number of outflows. The volume 
filling factor rises to a peak at a ~ 140° but is essentially constant from a — 110° — 180°. As indicated above, 
the volume per outflow is not conserved with varying opening angle and the volume of the enriched region 
actually decreases for increasingly anisotropic outflows. As a result one would naively expect the volume 
fllling factor to fall for decreasing a; however, since the number of galaxies increases the volume filling factor 
holds constant until a ^ 110° 

Figure [7] indicates that the number density of peaks hit decreases from 180° — 20°, but so does the 
volume filling factor (Fig. [5]). We can use both these statistics to determine the impact of varying alpha on 
the change in the number of peaks hit resulting from the path of least resistance of the outflows. We do this 
by taking the ratio of the number density of hits to the volume filling factor, shown in Figure [51 this is the 
number of hits per volume covered by outfiows. This indicates that as outfiows become more anisotropic 
they tend to avoid high-density structures and favor voids until a « 50° when this statistic turns around, 
indicating that most outfiows have crossed voids and have reached the next overdense structure. 



These statistics indicate that galactic outfiows undergo a transition over the range a — 100° — 50°, from 



- 23 - 




log [/)//?] 

Fig. 9. — The number of enriched grid points N'^ in the simulation vohime at z = 2 as a function of IGM 
overdensity for varying opening angle. In the top panel this is shown as a fraction of the number of systems 
at this density, Np. The statistic as a function of the total number of points, N, is shown in the middle 
panel. The bottom panel shows the number of enriched points below an overdensity threshold, N'p,^p by 
comparison with the isotropic case using $p given in equation (29). 

enriching their own high-density sources and the surrounding voids, to enriching neighboring high-density 
regions after crossing the voids. The top right panel indicates that in this range the radius of the largest 
outflow increases from around 0.3 — 0.4/i~^Mpc compared to the isotropic case of 0.23/i~^Mpc, confirming 
the plausibility of this explanation. 

4.2. Enrichment of Overdense and Underdense Systems 

We have established that the volume filling factor of enriched regions is dependent on the anisotropy 
of galactic outflows, but what is the nature of regions enriched? We investigate this by returning to the 

original power spectrum, and smoothing on a Jeans length scale to determine the baryonic density field 
produced with the same Gaussian random realization of structure. This density field is a result of purely 
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linear evolution from initial conditions with a Gaussian probability distribution function (PDF) . We map to 
a lognormal PDF in order to mimic a degree of non-linear behavior. This has been found to reproduce the 
PDF of the density in SPH simulations (|Bi fc Davidsen.il997i ) but is limited by an incomplete description of 
clustering. 

Having scanned through all 512"^ grid points in the mapped density field, we determined whether any of 
the galactic outflows reach them. Where this is the case we have flagged them as enriched and noted their 
overdensity. In Figure [5] we show three statistics derived from this approach for a range of different opening 
angles from 60° — 180°. We show the number of grid points enriched at a given overdensity, N'p, as a fraction 
of the number of grid points at that density, Np (top panel) ; this is effectively the probability of enriching a 
systems of a given density. This panel indicates that the effect of galactic outflows on overdense systems is 
dramatically reduced for increasingly anisotropic outflows. This is quite plausible in the context of the low 
number density of overdense systems. It is also clear that for decreasing opening angle the probability of 
enriching low-density systems increases. 

The middle panel shows the number of grid points enriched at a given overdensity as a fraction of the 
total number of grid points, N. In a plot of \og{p/p) these curves are Gaussian reflecting the underlying 
density PDF. The mean of the Gaussian, drifts to lower densities for lower opening angles. Even in the case 
of isotropic outflows the majority of metals are in underdense regions when considered by volume and as 
the outflows become anisotropic this effect becomes stronger as winds expand preferentially into low-density 
regions. 

The bottom panel is a cumulative version of the middle panel and so indicates the number of enriched 
grid points with an overdensity below a given threshold, N'^i^^. This statistic is shown by comparing to the 
isotropic case using 

'^P = K'<p/K'<PAS0- (29) 
This panel highlights the significant impact on the enrichment of underdense systems of anisotropic outflows. 
Anisotropic outflows can lead to an increase in the enriched volume of underdense systems of 10% (where 
a = 100° — 120°) and an increase of 40% in systems below p/p = 0.1 (where a = 80° — 100°) compared to 
isotropic outflows. The real volume filling factor can be determined and reaches a maximum of 11% — 12% 
for a w 140°. 



DISCUSSION 



There is a degree of uncertainty for the parameters /w, /dcpi /osc and they are of varying impor- 
tance. The star formation efficiency, has a well-documented and large uncertainty. It has a significant 
affect on the predicted volume filling factor as outlined in SFM. Our quoted values for the volume filling 
factor should be considered with this in mind. However, the dependence on opening angle will remain un- 
changed. The fraction of total supernova energy that contributes to the outflow, /w, and the mass escape 
fraction, /esc, both have associate d uncertainties. These uncertain ties have been theoretically constrained by 
SFM (and references therein) and Mori. Ferrara. &: Madaul 1 2002h . We are aware of no dir e ct observational 



constr a ints; although our r a tio of fesd f* is con sistent with observation bv iMartinI (jl999f ). iHeckman et al 
|2000l) . IPettini et al .1 |2000l) . iRupke et all |2005h and others. 



The value for the efficiency with which outflow metals are deposited upon collapsing density peaks, 
/dop, is poorly constrained and requires further analysis. This is not a concern here as the metal enrichment 
of halos does not have a significant impact on the subsequent numbers of galaxies and volume of the IGM 
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enriched. Only a handful of additional galaxies form before z = 2, with the inclusion of increased cooling 
by metal enrichment, that would otherwise not have formed (i.e. if /dop = 0). All density peaks that are 
enriched, and that will collapse to form galaxies in our simulations, will form earlier. This raises the distance 
traveled by outflows whose expansion is halted either by merger of the source galaxy or the end of the 
simulation nm. The volume enriched falls by less than 0.1%, when reducing the deposition efficiency from 
/dop = 0.9 to /dcp = 0. 

Where the outflows stall before their radii are larger than the pre-coUapse scale of their density peaks, 
we can reasonably assume that the outflows will fall back and galactic fountains will be formed. Since we 
allow all our outflows to join the Hubble flow when the internal pressure equals the external pressure, we 
have not attempted to simulate this effect. This will not change our results to a significant degree as this is 
the case in a small minority of our outflows. Even in those few cases, no stripping or metal enrichment of 
other peaks occur and they fill only a small fraction of the universe so making little impact on the volume 
filling factor statistics shown. 

In the calculation of external pressure we assume helium is singly-ionized, which is true for most of 
the redshift range where the outflows affect the IGM and will have little effect on the pressure. We also 
assume that the IGM is at the mean density and temperature and that the mean temperature is not redshift 
dependent. The density of the IGM is also implicitly assumed to be constant in the second term of equation 
(4) corresponding to the drag due to the sweeping up of the IGM (see TSE for more details). Most of the 
enriched volume is in underdense regions for all opening angles considered as shown in Figure 9 (discussed 
below) and the temperature of the IGM will be lower than that at mean density. Both these factors lower the 
external pressure, and so we can expect outflows to travel further after this correction. Outflows that only 
travel a short distance, such that they do not escape overdense regions, will travel even shorter distances as 
a result. A further complication arises since the outflows themselves will raise the temperature of the IGM. 
We will seek to investigate these issues in future work. 

We assume that the IGM is ionized in our calculation of the Compton drag on expanding outflows. 
This appears to b e reasonable since recent WMAP results indicate that the epoch of reionization is 2: ~ 11 



(|Page et al.ll2006f ) while outflows before z ~ 9 have a negligible impact on our results. 



The results from the previous section consistently indicate that we successfully describe anisotropic 
outflows that take a path of least resistance out of overdense regions, such as pancakes and filaments, and 
into voids. If sufficiently anisotropic (a < 100°) these outflows will also begin to strike neighboring overdense 
structures and further peaks. The number of peaks hit drops for anisotropic outflows and diminishes the 
ability of these outflows to strip peaks of their baryons. This reduces the capacity of ram pressure stripping 
at high- 2: to explain the abundance of dwarf galaxies in the Local Group. 

We find a value for the volume filling factor of 11% — 12% (for a « 140°), but this is sensitive to a 
number of factors such as the star formation efficiency and the degree of clustering of source galaxies. We 
will present a more thorough analysis of the value of this quantity in future papers in this series. We do, 
however, expect its dependence on opening angle of outflows to be rigorous. 

Despite this small volume filling factor, we find significant enrichment of underdense systems (below the 
mean density of the Universe), particularly when these outflows are anisotropic. The volume filling factor of 
enriched underdense systems is maximized by an opening angle of 100° — 120°, this is ~ 10% higher than that 
for isotropic outflows, while the volume of enriched space below p/ p = 0.1 is 40% higher. This may provide 
an explana tion for observation s (at high-z) of metal enrichm ent in systems at, or around, the mean density 



as seen in ISchave et al.l (|2003l ) and iPieri &: Haehnelti (|2004[) without the need to appeal to large volume 
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filling factors. These outflows also reach larger distances from the ir source galaxies and help explain why 
metal enrichment is seen far from observed galaxies at high redshift ( Fieri. Schave. fc Aguirrej[ 2006VSongailal 
2006.'). while still showing elevated metal enr ichment close to those galaxies ([Adelberger et al. .2003, .200. 



Fieri. Schave. fc: Aguirre 2006t Songaila 2006h . Since parts of the IGM are enriched solely based on whether 
they are located in the path of least res istance of an outflow , this may provide a fu rther source of scatter in 



the m etallicity of the IGM observed bv lSchave et al.l (|2003r ): ISimcoe et al.l (|2Q04I ): Fieri. Schave. &: Aguirre 
( 2006r) and as yet unexplained in simulations. In a subsequent paper in this series (jFieri. Grenon. & : MarteJ 
20071 ) we will investigate these issues by performing a direct comparison between these observations and 
synthetic QSO absorption spectra produced using our analytic description of anisotropic outflows. 

The strongest and weakest points of our method are the high dynamical range and the lack of gravita- 
tional dynamics, respectively. On the one hand, the ratio of the largest to smallest mass scale we consider is 
Afio/A/i « 50, 000. It would be very challenging for a numerical simulation to achieve such a large dynamical 
range in mass (that is, simulating galaxies and dwarf galaxies together), while having sufhcient resolution 
to properly simulate the outflows originating from the smallest galaxies. On the other hand, the treatment 
of large-scale structure and galaxy formation in the Monte Carlo method is quite simplistic. We combine 
a Gaussian random density field with a spherical collapse model for galaxy formation. In this approach, 
galaxies form at the comoving locations of density peaks and remain at these locations afterward. Even 
though we have a prescription for destruction of galaxies by mergers, the actual clustering of galaxies is not 
taken into account. If galaxies were allowed to cluster, collision and stripping by outflows would be more 
frequent, and also it would become more difficult to enrich low-density regions with metals. This lack of a 
correct description of dynamics means that the description of clustering of the IGM is also limited. 

The main limitation of this work is not the outflow model, but rather the Monte Carlo model used for 
describing galaxy formation. This Monte Carlo model has been used as a test bed for the outflow model 
i n order to perform an investigation of its importance and p otential influence. In two forthcoming papers 



( Martel. Grenon. fc Fieri 2007 : Fieri. Grenon. fc MarteJ l 120071 ). we will replace this Monte Carlo model by a 



more realistic numerical simulation of galaxy formation in a cosmological volume. 



6. SUMMARY AND CONCLUSION 

We have designed an analytical model for anisotropic galactic outflows based on the hypothesis that such 
outflows are bipolar and follow the path of least resistance through the environment of their source. In this 
analytical model we vary one parameter: the opening angle, a. We combined this model with an analytical 
Monte Carlo method for simulating galaxy formation, galaxy mergers, and supernova feedback. With this 
combined algorithm, we study the evolution of the galaxies and the IGM inside a comoving cosmological 
volume of size (12/i^^Mpc)^, from redshifts, z = 24 — 2, in a ACDM model. Our main results are the 
following: 

• Galaxy formation starts at redshift, z ~ 18. Since we neglect the formation and evolutionary times 
of massive stars, each newly formed galaxy immediately produces an outflow that lasts for a time, 
iburst ^ 50Myr. Such outflows can travel hundreds of kiloparsecs, and eventually collide with other 
objects. We neglect the effect of a collision with a well- formed galaxy (the cross-section is too small). 
When an outflow collides with a peak still in the process of collapsing, removal of the gas by ram 
pressure, preventing the formation of the galaxy, occurs about half of the time. When stripping does 
not occur, the protogalaxy is enriched in metals. This process occurs for all opening angles and the 
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proportion stripped or metal-enriched is essentially independent of opening angle. 

• When metal-enrichment of a peak occurs and this peak collapses to form a halo, the cooling time of 
that halo is reduced, and the galaxy forms earlier. However, this effect is rather small. In particular, 
we did not find that metal-enrichment could "bring to life" low-mass protogalaxies whose cooling time 
exceeds the age of the universe. 

• Anisotropic outflows channel matter preferentially into low-density regions, away from the cosmological 
structures (filaments or pancakes) in which the galaxies producing the outflows reside. Consequently, 
the number of halos encountered by outflows decreases with decreasing opening angle. This reduction 
in the number of hits results in a larger number of galaxies forming, since fewer halos are stripped by 
the ram pressure of outflows. 

• The volume filling factor of galactic outflows (that is, the volume fraction of the IGM occupied by 
outflows) holds constant and then decreases with opening angle. For angles a = 180° — 110° the 
constant filling factor is a result of the balance between an increase due to larger numbers of galaxies 
and a decrease due to a fall in volume per outfiow. At smaller angles, the volume of individual outflows 
drops significantly with a, and the total filling factor decreases since this term wins out. 

• The decrease in filling factor with decreasing angle is not sufficient to explain the decrease in number of 
hits. The ratio (number of hits) / (filling factor) decreases with decreasing angles down to a ^ 50°. This 
indicates that at these angles, the outflows are efficient at avoiding collisions with halos and channel 
matter preferentially into low-density region. Hence, if several halos reside in a common cosmological 
structure, an outflow produced by one of them will tend to avoid encountering the others. For angles 
a < 50°, we observe the opposite trend: outflows become more efficient in finding halos and hitting 
them. These narrow outfiows can travel across cosmological voids and hit halos located in unrelated 
structures, like the next filament or pancake. This effect is a continuation of a process begun at 
a ~ 100° where the mean distance travelled by outflows when they strip collapsing peaks of their 
baryons begin to increase as the first neighboring structures are hit. 

• The enrichment of the IGM with metals favors high-density systems since the sources of outflows are 
located in high-density regions. However, as the opening angle decreases, there is a dramatic reduction 
of the enrichment of such systems, combined with a dramatic increase in enrichment of low-density 
systems. Anisotropic outflows enrich around 10% larger a volume of the underdense Universe (and 
40% more of the Universe below p/p — 0.1) than isotropic outflows. 

This collection of results is a mere consequence of the fact that outflows follow the path of least resistance. 
This is an assumption in our model and is motivated by observations as well as high-resolution simulations. 

This work benefited from stimulating discussions with A. Ferrara, E. Scannapieco, J. Silk, M. Tegmark, 
R. Thacker, and S. D. White. All calculations were performed at the Laboratoire d'astrophysique numerique, 
Universite Laval. Figure [2] was produced at the Center for Computer Visualization, University of Texas, by 
Marcelo Alvarez. We thank the Canada Research Chair program and NSERC for support. CG is also 
supported by a Hubert Reeves fellowship. 
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A. NUMERICAL SOLUTION FOR THE OUTFLOW 

The equations governing the evolution of the outflow are 

n^H^R R^ ' 2 R^ ' ^ ' 

^ 27ri?3[i _cos(q!/2)] R ^ ' 

We solve these equations numerically, starting with the initial condition i? at i = 0, which we have 
chosen for conceptual simplicity. We might have chosen our initial radius at t = to be any value below 
the physical size of the star-forming region of the source galaxy; however, the results arc not substantially 
changed by such an adjustment. Also it is more plausible to expect that the starting point of the coherent 
outflows is closer to i? = than the radius of the star-forming region, if we assume that the star formation 
density will be centrally peaked. 

Most terms in equations (Al) and (A2) diverge at i = 0, making it impossible to obtain a numerical 
solution in this form. We must first perform a change of variables that will eliminate the divergences in 
these equations while retaining all the terms. To find the appropriate change of variables, we investigate 
the early-time behavior of the solution, by taking the limit i — > 0. In this limit, we have i? — > 0, -ff — > i/i, 
rib — > ^b,i^ — > fii, and L — s- Li. The gravity term, Hubble terms, and external pressure term are negligible 
(in the sense that they diverge slower than the other terms), and equation (IAl[l reduces to 



We can easily show that the solutions of the system of equations (jA2|) - (jA3[) are power laws, 

R = C•^^/^ (A4) 
_ 21fVH|C2 

^ ' 200^G ^ ' ^^^> 



where 



SOOGL, 



1/5 



C = <^ 77 , } and t->Q. (A6) 

' 23117b,,iJ2[l -cos(a/2)] ' ^ ' 



Using equations (A4) and (A5), we can find the proper change of variables. We introduce the following 
transformations (which are free from the above approximations), 

S = Rt^l^ , (A7) 
q = p^^/^ (A8) 
U ^ St. (A9) 

In the limit i ^ 0, the three functions S*, q, and U vary linearly with t. We stress that this leaves the 
behavior of the outflow expansion unchanged and will only lead to a well-behaved system of equations at 
the start of the first time step. We now eliminate the functions R and p in our original equations (Al) and 
(A2) using (A7)-(A9), and get 

19(7 Lt^ 5qU 



5t 27r[l - cos(a/2)]53 St 



(AlO) 
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9U US 

^2 



nhH^s 



5 



(All) 



(A12) 



These equations are completely equivalent to our original equations (Al) and (A2), but all divergences at 
t — have been eliminated. These equations can therefore be integrated numerically using a standard 
Runge-Kutta algorithm, with the initial conditions S = q = U = Oa,tt = 0. However, before doing so, it is 
preferable to rewrite the equations in dimensionless form. We define 



T = Hot, (A13) 

5* = HoS/C, (A14) 

U = HqU/C, (A15) 

q = Gq/HoC\ (A16) 

fH = H/Ho, (A17) 

h = L/L,, (A18) 

M = GM/C^Hy\ (A19) 



Equations (jA10p - (jA12[) reduce to their dimensionless form 

dq _ 19q 23inb,^fLfH,iT^ bqlJ 

~ IOOOttS'S 'J7' 

dS_ _ U_ 

dr T ' 

dU _ 9U US 87r(g-gcxt) 3t ( U 2S -V ^f^Sr 



dr 5t 25r ' Q^ff^S S \ t 5t ' j ' 2 

S^' 



(A20) 
(A21) 



(A22) 



where fH,i — Hi/Ho. Interestingly, the change of variable eliminates the explicit dependence upon the 
opening angle in equation (|A20p . 

The quantity /l appearing in equation (jA20|) depends on the luminosity Lcomp, which is given by 
equation We ehminate p and R in equation ^7^, using equations (jAT]) . (jAS]) . (|Al3p . (IAl4j) . and (IAl6p . 
then eliminate C using equation (jA6p . We get 

^comp ^ SOOtt'^ ath / kT^o \'^ z)^^^^ (A23) 
L, 2079 meHonb,^ffl^, \ He J ^ ' ^ ' 

The initial conditions arc S^q = U — OatT — 0. In the limit t ^ 0, where many terms take the form 
0/0, the derivatives reduce to dS/dr = dU /dr = 1, and dq/dr — 21i^b,ifH ^/SOOtt. We solve these equations 
numerically, using a fourth-order Runge-Kutta algorithm. 
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